Decorrelation of signals

ABSTRACT

A method of strongly decorrelating signals comprising processing input signals to determine delay and rotation parameters which implement at least one elementary paraunitary matrix. The parameters transform the input signals into output signals with improvement in a measure of strong decorrelation. The improvement in the measure of strong decorrelation is then assessed: if it is significant the latest output signals are, designated as input signals and the process is iterated. Iteration continues until the measure of strong decorrelation ceases to improve significantly, at which point the latest output signals are designated as signals decorrelated in a wide sense.

This invention relates to decorrelation of signals; more particularly itrelates to decorrelation of signals that have undergone convolutivemixing, and to a method, an apparatus and a computer program forimplementing decorrelation.

Decorrelation is known: it is a form of signal processing implemented bysoftware running on a computer system and accepting data from sensorsafter conversion from analogue to digital form. Two assumptions arenormally made in decorrelation, stationarity and linearity, and theseassumptions are also made in connection with the present invention.Stationarity means that signals and channels in which they mix do notchange over a time interval during which mixed signals are sampled.Linearity means that mixtures of signals received by sensors are linearcombinations of these signals. More complicated combinations featuringsignal products and squares and higher order powers of signals are notconsidered.

The aim of decorrelation is to remove similarities between inputsignals, thereby eliminating redundancies. In some circumstances thismay lead to recovering the signals as they were prior to mixing, butmore usually it is used as a pre-processing step in an unmixingalgorithm, or as part of another algorithm such as a data compressionalgorithm. A common usage is to separate the signals which carrymeaningful data from those that arise due entirely to noise in theenvironment and in a receiver system which detects the signals: suchseparation may be an end in itself, or it may be used as part of a datacompression scheme.

Decorrelation is a process of removing correlations or similaritiesbetween signal pairs in a set of signals: correlation of a pair ofsignals is itself defined mathematically as an integral of the signalpair's product over time. It is important that as little information aspossible (preferably none) is destroyed in decorrelation process, aswould occur for example if decorrelation led to one of the signals beingzeroed. Use of energy preserving transforms for decorrelation avoidsinformation being destroyed and leads to output signals with otheruseful properties.

If an energy preserving transform is used, it is usually required torecover decorrelated signals in decreasing order of power. This allowsdata compression techniques and signal/noise subspace separationtechniques to be implemented very simply.

Algorithms for decorrelating signals that have been instantaneouslymixed are known: they may use either transforms that are energypreserving or those that are not.

Unfortunately, an algorithm which is adequate for decorrelating signalsthat have been instantaneously mixed cannot cope with more difficultproblems. These more difficult problems occur when an output signal froma sensor must be expressed mathematically as a convolution, i.e. acombination of a series of replicas of a signal relatively delayed withrespect to one another. It is therefore referred to as the “convolutivemixing” problem.

The approach used in instantaneous algorithms has been extended to theconvolutive mixing situation: from this approach it has been inferredthat convolutively mixed signals can be decorrelated by the use of anenergy preserving technique. This algorithm would accommodate timedelays involved in mixing and decorrelation. The criterion to be imposedhas however changed: i.e. pointwise decorrelation (as defined later) isno longer sufficient, instead a wider condition involving decorrelationwith time delays should be imposed. This will be referred to as widesense decorrelation, and imposing this property as imposing strongdecorrelation.

Instead of the unitary matrix employed in instantaneous algorithms, analgorithm for imposing strong decorrelation by an energy preservingfilter employs a paraunitary matrix. As will be described later in moredetail, a paraunitary matrix is one which gives the identity matrix whenmultiplied by its paraconjugate matrix—a polynomial equivalent of aHermitian conjugate matrix. A possible approach for the strongdecorrelation problem is therefore to search for a paraunitary matrixthat maximises a measure of the wide sense decorrelation.

In this case, it may be possible to improve on the mere return ofsignals in decreasing order of power. It may be possible to provide forreturn of signals in decreasing order of power at each frequency in aset of frequencies. Thus a first signal has more power than othersignals at all frequencies and a second signal has more power than allother signals apart from the first at all frequencies etc. This propertyis called spectral majorisation. Spectral majorisation is useful if thealgorithm is to be used for data compression or for separating signaland noise subspaces.

In “Theory of Optimal Orthonormal Filter Banks” ICASSP 1996, P. P.Vaidyanathan discloses a use for filtering techniques where both strongdecorrelation and spectral majorisation are necessary, neither beingsufficient without the other.

The simplest prior art method of imposing strong decorrelation involvesusing a multichannel whitening lattice filter. See for example S.Haykin, “Adaptive Filter Theory”, Prentice Hall, 1991, although thereare many references for these techniques, based on slightly differentalgorithms designed with different convergence and stability propertiesin mind. These techniques were not designed to impose strongdecorrelation directly: instead they aim to recover an innovationssequence of the mixed signals. The components of the innovationssequence so produced are strongly decorrelated as a consequence of whatthey are. However these techniques are not energy preserving, and so arenot suitable for use in many scenarios.

In “Principal Component Filter Banks for Optimal MultiresolutionAnalysis”, IEEE transactions on Signal Processing Vol 43 Aug. 1995, M.K. Tsatsanis and G. B. Giannakis discuss how to find an optimal energypreserving and spectral majorising filter. This reference uses afrequency domain technique to show how an optimal filter can be found ifinfinite information is held about the signals and no constraint isplaced upon the order of the filter. However it does not show either howthis can be transferred to a data dependent algorithm, or how the filtercan be constrained to a sensible order without losing optimality.However its calculations provide absolute performance bounds, and formedthe basis of several of algorithms discussed below.

In “Multirate Systems and Filter Banks”, Prentice Hall: SignalProcessing Series, 1993, P. P. Vaidyanathan discloses parameterisationof paraunitary matrices in a stage-by-stage decomposition of aparaunitary matrix in z⁻¹: here z⁻¹ is an operator implementing a delay.Vaidyanathan shows that a product matrix built up from a series of pairsof paraunitary matrix blocks is paraunitary: here, in a pair of blocks,one block represents a delay and the other a 2 by 2 unitary matriximplementing a Givens rotation (see U.S. Pat. No. 4,727,503).Vaidyanathan also shows that a para unitary matrix of degree N, definedlater, is a product of N+1 rotations and N one-channel delay operatorsall implementing the same unit delay.

In “Attainable Error Bounds in Multirate Adaptive Lossless FIR Filters”,ICASSP 1995, P. A. Regalia and D Huang propose two differenteigenvalue-based algorithms for strong decorrelation. Their algorithmsare recursive, limited to two input signals and based upon a fixeddegree lossless filter with parameters disclosed by Vaidyanathan. It isalso necessary to access data from various points internal to theVaidyanathan decomposition for use in adjusting the parameters: thisadds an extra layer of processing complexity.

In “Rational Subspace Estimation using Adaptive Lossless Filter” IEEEtransactions on Signal Processing Vol 40 October 1992 P. A. Regalia andP. Loubatton disclose a different algorithm for attempting to imposeboth strong decorrelation and spectral majorisation. This algorithm isapplicable to more than two signals, but it is based upon a deflationaryapproach where an optimal single signal output is first calculated, andthen augmented with other outputs until one of the latter becomesinsignificant. Again a fixed degree and fixed order lossless filter isparameterised as disclosed by Vaidyanathan, and the parameters areadjusted according to a stochastic algorithm, which requires informationon the data flow inside the Vaidyanathan decomposition.

A different algorithm is proposed in “Design of Signal-Adapted FIRParaunitary Filter Banks” ICASSP 1996 by P. Moulin, M. Anitescu, K. O.Kortanek and F. Potra. This uses a slightly different version of theVaidyanathan parameterisation with a fixed degree. From the effects of aparameterised filter on one single output channel, this algorithm showshow to convert a problem of maximising power in one channel into alinear semi-infinite programming problem, which has known methods forsolution. In effect the algorithm has an objective of achieving bothstrong decorrelation and spectral majorisation; it has a cost functionwhich is very similar to a known N₁ ⁰ cost function to be described inmore detail later. However it forces the paraunitary filter to be of afixed degree, and relies upon a link between parameters and costfunction being computationally tractable. For different cost functionsthis would not be the case.

Gradient descent methods are known which aim to adjust all parameters ofa paraunitary matrix or an unmixing filter simultaneously: these methodshave another difficulty, i.e. they have parameters linked to any usefulmeasure of independence in a very complex way which does not factoriseeasily. This means adjusting all parameters at once: it leads to veryslow algorithm convergence, which is why none of the algorithmspreviously mentioned use gradient descent. Gradient descent tends to beavoided in the prior art due to problems with speed and accuracy ofconvergence.

It is an object of the invention to provide an alternative method ofstrongly decorrelating signals.

The present invention provides a method for strong decorrelation ofsignals which are input, characterised in that it includes the stepsof:—

-   a) processing the input signals to determine delay and rotation    parameters which implement at least one elementary paraunitary    matrix and transform the input signals into output signals to obtain    improvement in a measure of strong decorrelation;-   b) assessing the improvement in the measure of strong decorrelation,    and if it is significant designating the output signals as input    signals and iterating step a) and this step b);-   c) if the improvement is not significant designating the output    signals as signals decorrelated in a wide sense.

The invention provides a technique for wide sense decorrelation ofsignals which have been convolutively mixed, even if mixed using unknownfilters. It is suitable, for example, for signals which have undergonemultipath reflection, or acoustic signals in a reverberant environment.It is believed that the method of the invention is more reliable andrequires less calculation than any alternative currently known. Byvirtue of use of an elementary paraunitary matrix, decorrelation iscarried out using an energy preserving transform, and examples of theinvention may (not necessarily but if required) impose spectralmajorisation on the signals as well as imposing wide sensedecorrelation.

To transform the input signals, the method of the invention maydetermine in step a) delay and rotation parameters which characterise asingle elementary paraunitary matrix. It may include producing aparaunitary matrix by cumulatively multiplying successive elementaryparaunitary matrices produced by iterating step a).

The range of signal delay parameters may be a set of discrete delayvectors, and the delay and rotation parameters may be determined bygenerating a respective version of the input signals delayed by eachdelay vector in the set, and for each version finding rotationparameters which at least approach producing maximisation of the outputsignals' strong decorrelation. The rotation parameters which at leastapproach producing maximisation of the output signals' strongdecorrelation may be determined using an algorithm for instantaneousdecorrelation. The at least one elementary paraunitary matrix mayimplement at least one leading delay, rotation and terminal delay.

The method may involve n input signals where n is an integer greaterthan 2, the range of signal delay parameters being a set of n-elementdelay vectors and the range of signal rotation parameters being a set ofn(n−1)/2 angle parameters. Step a) may comprise determining delay androtation parameters which implement at least one elementary paraunitarymatrix providing for rotation of a pair of input signals and relativedelay of the or as the case may be each other input signal. The n inputsignals may be associated with respective channels, step a) havingn(n−1)/2 successive stages each associated with at least one respectiveelementary paraunitary matrix and each providing for rotation of signalsassociated with a respective pair of channels and provision of relativedelay associated with the or as the case may be each other channel, thefirst stage being arranged to process the input signals and the or asthe case may be each subsequent stage being arranged to receive signalsprocessed in the respective preceding stage.

The method of the invention may involve a set of n input signals where nis an integer greater than 2, characterised in that it comprises:

-   a) producing n(n−1)/2 replicas of the set of input signals,-   b) in each replica selecting a respective signal pair differing to    that selected in other replicas, and-   c) the step of processing the input signals to determine delay and    rotation parameters being carried out for each replica and    comprising    -   i) determining delay and rotation parameters which implement at        least one elementary paraunitary matrix providing for rotation        of the respective selected signal pair only, and    -   ii) determining which replica when transformed by the associated        at least one elementary paraunitary matrix gives rise to        transformed signals corresponding to improvement in a measure of        strong decorrelation by at least a major part of a maximum        extent obtainable over the replicas and designating these        transformed signals as output signals.

In another aspect, the present invention provides computer apparatus forstrong decorrelation of signals, the apparatus being programmed forreception of input signals, characterised in that the apparatus is alsoprogrammed:—

-   a) to process the input signals to determine delay and rotation    parameters which implement at least one elementary paraunitary    matrix and transform the input signals into output signals to obtain    improvement in a measure of strong decorrelation;-   b) to assess the improvement in the measure of strong decorrelation,    and if it is significant designate the output signals as input    signals and iterating a) and b);-   c) if the improvement is not significant to designate the output    signals as signals decorrelated in a wide sense.

In a further aspect, the present invention provides a computer programmefor implementing strong decorrelation of signals input to computerapparatus, characterised in that it has instructions for controllingcomputer apparatus:—

-   d) to process the input signals to determine delay and rotation    parameters which implement at least one elementary paraunitary    matrix and transform the input signals into output signals to obtain    improvement in a measure of strong decorrelation;-   e) to assess the improvement in the measure of strong decorrelation,    and if it is significant designate the output signals as input    signals and iterating a) and b);-   f) if the improvement is not significant to designate the output    signals as signals decorrelated in a wide sense.

The apparatus and computer programme aspects of the invention mayinclude features equivalent to those mentioned above in connection withthe method of the invention.

In order that the invention might be more fully understood, anembodiment thereof will now be described, by way of example only, withreference to the accompanying drawings, in which:—

FIG. 1 illustrates decomposition of a paraunitary matrix into degree-oneparaunitary matrices in a prior art convolutive unmixing process;

FIG. 2 is a schematic block diagram of decomposition of a paraunitarymatrix into elementary paraunitary matrices in a convolutive wide sensedecorrelation process of the invention;

FIG. 3 is a schematic block diagram illustrating production of aparaunitary matrix from elementary paraunitary matrices in a process ofthe invention;

FIG. 4 illustrates production of elementary paraunitary matrices in aprocess of the invention;

FIG. 5 illustrates decomposition of a paraunitary matrix into elementaryparaunitary matrices in a prior art wide sense decorrelation processinvolving more than two signal channels;

FIGS. 6 and 7 schematically illustrate alternative wide sensedecorrelation processes of the invention involving more than two signalchannels;

FIG. 8 illustrates implementation of an alternative form of elementaryparaunitary matrix involving a leading delay, a rotation and a terminaldelay; and

FIG. 9 illustrates use of the FIG. 8 elementary paraunitary matrix in anequivalent of the FIG. 4 process.

The present invention is arrived at by modifying and extending the priorart of both instantaneously decorrelating signals and of stronglydecorrelating signals: this prior art will first be discussed in moredetail to enable the invention to be better understood.

Wide sense decorrelation is based entirely upon second order propertiesof data. The equality at (1) below holds if and only if theinstantaneous decorrelation property holds for the pair of signals x₁,x₂.E(x ₁ x ₂)=E(x ₁)E(x ₂)  (1)

Here E is the expectation operator, and a product x₁x₂ represents thesignal obtained by multiplying each sample of a signal x₁ by arespective sample of a signal x₂ having the same sample time: thisdefines pointwise decorrelation of a pair of signals.

Expectations can be estimated from sampled signals to determine whetheror not Equation (1) holds for sampled signals. If Equation (1) holds forall signal pairs then the signals are referred to as being decorrelated.

Mathematically, instantaneous mixing of two signals from sources (e.g.loudspeakers) and their receipt by two receivers can be represented asfollows:—x ₁(t)=h ₁₁ s ₁(t)+h ₁₂ s ₂(t)x ₂(t)=h ₂₁ s ₁(t)+h ₂₂ s ₂(t)  (9)

Here loudspeaker outputs (original signals) reaching the two receiversat time (t) are represented by s₁(t) and s₂(t): t is in units of aninterval between successive digital signal samples (clock cycles), andhas integer values 1, 2 etc. The path from each loudspeaker to eachreceiver and the receiver's reception characteristics give rise to aconstant channel modulation. In the instantaneous mixing problem, thismodulation is represented by multiplier terms h_(ij), where i=1 or 2indicates first or second receiver and j=1 or 2 indicates first orsecond loudspeaker respectively. First and second loudspeaker signals attime (t) are represented by x₁(t) and x₂(t). Each receiver signal is thesum of two loudspeaker signals, each multiplied by a respectivecoefficient. These equations can be represented more compactly inmatrix-vector format as:—x(t)=Hs(t)  (3)where s is a 2 by 1 element vector of two loudspeaker signals, H is a 2by 2 element mixing matrix and x is a 2 by 1 element vector of tworeceived signals, one at each receiver. H is referred to as the mixingmatrix because it has matrix elements (coefficients) that operate on theloudspeaker signals and mix them to form the received signals. The sizesof s, x, and H will increase if more than two signals or more than twosensors are involved. As previously mentioned, algorithms for solvingthe instantaneous decorrelation problem produced by Equation (3) existin the prior art.

Mathematically, a correlation matrix R of a set of signals expressed asa vector x can be defined as:—R=E((x−E(x))(x−E(x))^(H))  (4)where E(x) is the expectation of x, that is the mean of x over time andthe superscript index H means a Hermitian conjugate. In Equation (4)E(x) is subtracted from x, which is the first step in all processingthat follows. However, to simplify notation it is usually assumed thatsignals have zero mean values, and this assumption will be used in whatfollows. This means that Equation (4) simplifies to R=E(xx^(H)).However, it is not necessary for signals to have zero mean values inorder for the invention to be effective.

To impose decorrelation upon mixed signals, it is necessary to force thecorrelation matrix R to have off-diagonal coefficients which are equalto zero. This is the instantaneous decorrelation problem. It can bedone, for example, by ensuring that a correct amount of received signal1 is subtracted from received signal 2 to ensure that the resultingoutput is decorrelated from received signal 1. This is repeated for allother signals, producing a set of outputs decorrelated with receivedsignal 1, which is then defined as an output 1. The remaining outputsare then decorrelated with each other by repeating this procedure forthem. It is equivalent to obtaining decorrelated signals x′ bypre-multiplying signals x by a lower triangular matrix L, i.e. a matrixhaving all above diagonal coefficients equal to zero. Instantaneouslydecorrelated signals are designated x′ and have a correlation matrix R′defined by:—R′=E(x′x′ ^(H))=E(Lxx ^(H) L ^(H))=D  (5)

Here the new correlation matrix is shown as D, which is diagonal andpositive semi-definite.

In order that instantaneous decorrelation algorithms be capable ofimposing second order independence, a relative propagation delay betweena time sample of a signal from a signal source reaching two separatesensors must representable as a simple phase shift applied to a signalreceived by one of the sensors. For N sensors this becomes N−1 phaseshifts. If this condition does not hold, channels from signal sources tosensors may be convolutive and in consequence instantaneousdecorrelation algorithms will fail to impose second order independence.

One-channel convolution and its analysis using z-transforms will now beconsidered. In a single signal, single receiver arrangement, theloudspeaker signal consists of a series of values indexed by therespective times t at which they were obtained and denoted by s(t). Theloudspeaker signal passes through a channel and is received by thereceiver that detects another series of values, x(t). The link betweens(t) and x(t) may not be a simple multiplication as in instantaneousmixing. Instead x(t) may consist of a linear sum of several of earliervalues of s(t). This is known as a convolution, and is shown by:—$\begin{matrix}{{x(t)} = {{h \otimes {s(t)}} = {\sum\limits_{k = 0}^{p}{{h(k)}{s\left( {t - k} \right)}}}}} & (6)\end{matrix}$where {circle around (×)} is referred to as a “convolution operator”.Generation of a linear sum of values of s(t) may occur because signalspass through an irregular channel which spreads them out. It may alsooccur from multi-path effects, i.e. when there is more than one path asignal can take to reach a receiver and the paths have differingpropagation times. Instantaneous mixing is a simple case of convolutivemixing.

Convolution is more complex and involves more algebra than simplemultiplication by a coefficient. However it is possible to reduce it toa multiplication by the use of z-transforms as follows. The z-transformsof s, x and h are:—s(z)=s(0)+s(1)z ⁻¹ +s(2)z ⁻² +s(3)z ⁻³+ . . .x(z)=x(0)+x(1)z ⁻¹ +x(2)z ⁻² +x(3)z ⁻³+ . . .h(z)=h(0)+h(1)z ⁻¹ +h(2)z ⁻² + . . . +h(p)z ^(−p)  (7)

Here z⁻¹ is referred to as a delay operator as mentioned above. Thez-transform of a signal expressed as a time series of sample values is apolynomial in z⁻¹. Given the form of the z-transform, its polynomialcoefficients enable recovery of the original signal. When in theirz-transform forms, the convolution of s by h as shown in Equation (6)becomes a multiplication of two polynomials. Thus:—x(z)=h(z)s(z)  (8)

Convolutive mixing may be addressed by combining instantaneous mixing asper Equation (1) with the more complicated approach of Equations (6) to(8). This deals with a plurality of sources and a plurality of receiverswith intervening channels which are convolutive. Thus for thetwo-signal, two-receiver case the following expressions for receivedsignals can be written:—x ₁(t)=h ₁₁ {circle around (×)}s ₁(t)+h ₁₂ {circle around (×)}s ₂(t)x ₂(t)=h ₂₁ {circle around (×)}s ₁(t)+h ₂₂ {circle around (×)}s₂(t)  (9)

In the matrix-vector notation of Equation (3), and using the z-transformnotation of Equations (7) and (8), Equation (9) can be written simplyas:—x(z)=H(z)s(z)  (10)

Here s(z) is a 2 by 1 coefficient vector formed from z-transforms ofloudspeaker outputs. H(z) is a 2 by 2 coefficient polynomial mixingmatrix and x(z) is the 2 by 1 vector formed from z-transforms ofreceived signals. Again if more than two loudspeakers and/or receiversare present, the sizes of s(z), x(z) and H(z) increase accordingly.

H(z) is called a polynomial mixing matrix because it represents mixingof signals: its coefficients are polynomials in z⁻¹. It may also beconsidered to be a polynomial in z⁻¹ with coefficients which arematrices. It is found one coefficient at a time by taking z-transformsof individual loudspeaker/receiver channels. For future convenience, theorder and degree of a polynomial matrix will now be defined. The orderof a polynomial matrix is the greatest power to which z⁻¹ is raised inthe matrix. The degree of a polynomial matrix is the smallest number ofdelays necessary to implement it as a filter. It is always at least thesame size as the order but may be greater.

The instantaneous decorrelation problem only deals with removingcorrelation between the signals with no time delays applied to them. Inorder to impose wide sense decorrelation the measure of decorrelationhas to operate across time delays. That is, to be strongly decorrelated,the correlation between two signals has to be zero when any time delayis applied to one of the signals. This can be written mathematicallyas:—R(d)=E(x(t)x(t−d ^(H)))=D(d)∀dε( . . . −2, −1, 0, 1, 2, . . . )  (11)where ∀dε( . . . −2, −1, 0, 1, 2, . . . ) means all values of delay d(in units of a time interval between successive signal samples) in theset of all positive and negative whole numbers, and D(d) is a diagonalmatrix (all off-diagonal coefficients zero) for all values of d. It ispossible to put Equation (11) into z-transform form as follows:—R(z)=E(x(z){tilde over (x)}(z))=D(z)  (12)where D(z) is a diagonal polynomial matrix. Equation (12) introduces theconcept of a paraconjugate {tilde over (x)}(z) of a polynomial matrixx(z) denoted by a tilde above the matrix symbol {tilde over (x)}. Itmeans combined operations of transposition and conjugation (as in theHermitian operator) and the substitution of 1/z* for z. Paraconjugationis an extension of the Hermitian operator to polynomial matrices.

To achieve wide sense decorrelation it is necessary to force alloff-diagonal terms of a polynomial correlation matrix to become the zeropolynomial, i.e. the polynomial in z⁻¹ having all its coefficients equalto zero. No constraints are put upon the polynomial functions of thematrix on its diagonal. Thus:— $\begin{matrix}{{R^{\prime}(z)} = \begin{bmatrix}{\sigma_{1}(z)} & 0 \\0 & {\sigma_{2}(z)}\end{bmatrix}} & (13)\end{matrix}$

The unit polynomial is the polynomial that has zero coefficients for allpowers of z⁻¹ except for z⁰ which has unity as its coefficient. In theprior art a multichannel lattice filter, see e.g. S. Haykin, “AdaptiveFilter Theory”, Prentice Hall, 1991, imposes strong decorrelation whileensuring that the diagonal terms of the polynomial correlation matrixare the same by forcing them all to be the unit polynomial. This whitensthe signals.

The polynomial correlation matrix of the signals x(z) is given byEquation (12) above. A multichannel lattice filter transforms these toobtain a matrix of strongly decorrelated and whitened signals x′(z).These signals have the property that their correlation matrix R′(z) isthe identity polynomial matrix I, symbolically:—R′(z)=E(x′(z){tilde over (x)}′(z))=I  (14)

The fact that this method imposes the whitening condition upon thesignals means that it carries out more calculation than is necessary andhas not necessarily used an energy preserving transform. It can alsolead to difficulties when original signals driving the mixing were notwhite, or when mixing nulls some of their frequency components, as itwhitens the outputs and thus amplifies the noise in those frequenciesthat were nulled.

Instantaneous algorithms that use energy preserving transforms achievedecorrelation by the use of unitary matrices. These are energypreserving, and are defined as those matrices that when multiplied bytheir Hermitian conjugate give the identity. The extension of theHermitian operator to the space of polynomial matrices is the operationof paraconjugation. This leads to the definition of paraunitary matricesas those polynomial matrices that give the identity matrix whenmultiplied by their paraconjugate. Symbolically H(z) is a paraunitarymatrix if and only if H(z){tilde over (H)}(z)={tilde over (H)}(z)H(z)=I.

Paraunitary matrices are energy preserving and, more strongly, areenergy preserving at all frequencies. This means that if decorrelationis achieved by the use of paraunitary matrices, the problems ofwhitening mentioned above are avoided.

As previously mentioned, the Vaidyanathan reference disclosesparameterisation of the space of all paraunitary matrices. This gives astage-by-stage decomposition of a paraunitary matrix in z⁻¹. Theparameterisation is illustrated in FIG. 1 for a two-channel case. Itcomprises a series of N+1 rotation blocks Q₀, Q₁, . . . Q_(N), adjacentpairs of these blocks being separated by individual delay blocks Λ(z)which are all alike. Upper and lower signal channels 20 and 22 passthrough all the blocks. Amplifiers 24 indicate channel input signalscaling factors α, of which the upper channel factor is positive and thelower channel factor is positive or negative as indicated by “+−” before“α”.

Each of the rotation blocks Q₀ to Q_(N) implements a respective 2 by 2unitary matrix which itself implements a Givens rotation (see U.S. Pat.No. 4,727,503): a Givens rotation is a rotation parameterised by asingle rotation angle θ. In the first rotation block Q₀, a signal in thelower channel 22 is multiplied by Givens rotation parameters s₀ and c₀;a signal in the upper channel 20 is multiplied by Givens rotationparameters −s₀ and c₀; s₀ and c₀ are respectively the sine and cosine ofthe rotation angle θ implemented by the block Q₀. Each c₀ product issummed with the s₀ or −s₀ product from the other channel to provide asum which passes along the respective channel 20 or 22 to the next blockΛ(z). This next block Λ(z) delays the signal in the lower channel 22with respect to that in the upper channel 20. This procedure is thenrepeated in subsequent rotation block Q₁ to Q_(N), with equalintervening delays Λ(z) in the lower channel 22.

The rotation blocks Q₀ to Q_(N) and delays Λ(z) are themselvesparaunitary, and so their product is a matrix which is also paraunitary.Vaidyanathan proved that any finite degree N paraunitary matrix can beexpressed in this form. Thus apart from a scalar factor, any finitedegree paraunitary matrix can be expressed as follows:—H _(N)(z)=Q _(N) . . . Λ(z)Q ₁Λ(z)Q ₀  (15)where N is the degree of the paraunitary matrix. Thus a paraunitarymatrix H_(N)(z) of degree N is the product of N+1 Givens rotations and None-channel delay operators of equal delays, as shown by theconcatenated stages in FIG. 1.

However, there is a difficulty with using Vaidyanathan parameterisation,a simple example of which is given here: consider signals which form anindependent identically distributed (iid) time series. Assume thesesignals have been mixed by a paraunitary mixing matrix. Let theparaunitary mixing matrix have degree 1 and be formed, usingVaidyanathan parameterisation, as:— $\begin{matrix}\begin{matrix}{{H_{N}(z)} = {Q_{1}{\Lambda(z)}Q_{0}}} \\{= {{\begin{bmatrix}1 & 0 \\0 & 1\end{bmatrix}\begin{bmatrix}1 & 0 \\0 & z^{- 1}\end{bmatrix}}\begin{bmatrix}{\cos\left( \theta_{0} \right)} & {\sin\left( \theta_{0} \right)} \\{- {\sin\left( \theta_{0} \right)}} & {\cos\left( \theta_{0} \right)}\end{bmatrix}}}\end{matrix} & (16)\end{matrix}$i.e. the last of the two Givens rotations (within the left hand pair ofsquare brackets) is in fact the identity. It has been discovered inaccordance with the invention that attempts to impose strongdecorrelation on signals mixed by H_(N)(z) should begin by undoing thedelay, not by applying a rotation. However, in applying a single blocktype approach, the first step is to try to find an appropriate rotationto apply, even though one is not needed. This superfluous rotation isvery difficult for later blocks to undo, and its application canintroduce more correlation between the signals, which is precisely theopposite of the required decorrelation.

The present invention is intended to overcome the difficulties indicatedabove, and to avoid superfluous calculations. It has similarities withVaidyanathan's parameterisation but unlike the latter it employsvariable delays not necessarily confined to one channel and notnecessarily equal; moreover it does not force a rotation to be appliedwhen there is little or nothing to gain, and allows arbitrarily sizedparaunitary matrices to be constructed.

FIG. 2 illustrates the construction of a paraunitary matrix H_(N)(z) asa product of N+1 elementary paraunitary matrices in accordance with theinvention. An elementary paraunitary matrix is an expression coined forthe purposes of the present invention: it is defined as a polynomialmatrix which represents applying a set of delays to the signals,followed by a single unitary transformation (rotation); as will bedescribed later in more detail, this definition may optionally includethe rotation being followed by a further set of delays (omitted in thepresent example). In FIG. 2 and as in FIG. 1, for the sake ofsimplicity, only the two-channel case is shown as indicated by upper andlower channels 30 and 31. Each of a set of dotted line blocks V_(i)(z)(i=0, 1, to N) represents a respective single elementary paraunitarymatrix that implements a respective delay d_(i) followed by a Givensrotation. N+1 of the blocks V_(i)(z) are used to make up the paraunitarymatrix H_(N)(z).

Each of the blocks V₀(z) to V_(N)(z) implements a respective 2 by 2elementary paraunitary matrix which itself implements a Givens rotation.In the first block V₀(z), signals in the upper and lower channels 30 and31 pass into a matrix of four ganged single pole double throw switches32 controlling routing or otherwise via a delay cell Z^(d0). Theswitches 32 are all either in the UP position or in the DOWN position:when DOWN, upper channel signals are delayed at Z^(d0) but lower channelsignals are not, and when UP the reverse is the case. This enableseither channel 30 or 31 to be delayed relative to the other.

On leaving the switches 32, an upper channel signal is multiplied byGivens rotation parameters −s₀ and c₀ at amplifiers 33 and 34respectively. Similarly, on leaving the switches 32, a lower channelsignal is multiplied by Givens rotation parameters s₀ and c₀ atamplifiers 35 and 36 respectively. Here s₀ and c₀ are respectively thesine and cosine of a rotation angle θ₀ implemented by the block V₀(z).Each c₀ product is summed with the so or −s₀ product involving thesignal from the other channel: this provides sums, which pass alongrespective channels 30 and 31 to the next block V₁(z). This delay/rotateprocedure is then repeated in subsequent blocks V₁(z) to V_(N)(z). Theselater blocks will not be described further because they operate in thesame way as the first block V₀(z), except their delays and rotationswill usually be different. Implementation of this procedure to findappropriate delays and rotations will be described later, following atheoretical discussion.

A full paraunitary matrix parameterisation in accordance with theinvention in terms of elementary paraunitary matrices can be expressedas follows:—H _(N)(z)=V _(d) _(N) _(,θ) _(N) (z) . . . V _(d) ₁ _(,θ) ₁ (z)V _(d) ₀_(,θ) ₀ (z)  (17)

The progress of imposing wide sense decorrelation upon signals inaccordance with the invention uses a cost function that is calculatedafter every application of an elementary paraunitary matrix to thesignals to determine whether or not an improvement has been achieved.The cost function is a measure of the signals' strong decorrelation, anexample of which is a sum squared of second order measures ofautocorrelation of the signals, and is at a maximum when allcross-correlation terms are at a minimum. Thus an increase in this costfunction corresponds to a reduction in correlation between the signals.The cost function will be referred to as the Decorrelation Cost Function(hereinafter “DCF”).

The process of the invention is to find an appropriate paraunitarymatrix that imposes wide sense decorrelation on the signals. This matrixis not parameterised as disclosed by Vaidyanathan, but it is insteadparameterised in terms of elementary paraunitary matrices (see Equation(17)) calculated sequentially. Once derived, each elementary paraunitarymatrix is applied to the signals used to calculate it: this transformsthe signals into a form suitable for calculating the next such matrix.

To calculate an elementary paraunitary matrix a set of possible delaysis considered. A latest elementary paraunitary matrix is chosen whichmaximises the DCF calculated for signals transformed by application ofall such matrices calculated up and including the latest. To identifythe latest elementary paraunitary matrix its delay and rotationparameters are required. There is a countable infinity of differentpossible delays, but an uncountable number of different rotationparameters. Realistically the countable infinity of possible delays willbe reduced to a finite number of different delays to be considered,either by constraining them or by making the assumption of circulantstatistics and hence constraining them to be less than the data length(number of samples of each signal used in decorrelation).

Once a trial value of delay has been applied to the signals, rotationparameters s₀ and c₀ are chosen. These parameters could be chosen byusing a technique which explicitly maximises the DCF of output signals.

It is therefore possible to look at a range of trial delays and for eachdelay to attach the rotation that maximised the DCF after theapplication of the delay and the rotation. This produces a set ofpossible elementary paraunitary matrices, each associated with arespective DCF found by applying it to signals and finding the DCF ofthe resulting output. Under certain assumptions it may be possible toavoid applying the elementary paraunitary matrices to the signals, andinstead to calculate the DCF of the resulting output directly. This hascomputational advantages, but is only possible in some cases.

The elementary paraunitary matrix that is selected is that having thebest DCF output. This stage is repeated until there is no gain in DCF tobe found, or the gain in DCF is regarded as being too small to warrantthe extra complexity in the strongly decorrelating matrix.

Thus the invention can be seen as a method of imposing wide sensedecorrelation by sequentially applying a series of elementaryparaunitary matrices each chosen to maximise an DCF measure in an outputresulting from its use. Each elementary paraunitary matrix contains onerotation. Thus the process of the invention implements an algorithmreferred to as the Sequential Best Rotation (SBR) algorithm. Thisalgorithm comprises decomposition of the paraunitary matrix intoelementary paraunitary matrices, and the use of a measure which can bemeaningfully calculated at each stage of the algorithm not just at theend. The combination of these two steps separates the algorithm into aseries of like stages that can be carried out in sequence.

Strong decorrelation of signals is equivalent to a process of generatinga strongly decorrelating matrix H(z). This process will now bedescribed. It is not essential to generate H(z), but it permits theprocess to be tested: i.e. mixed signals which have undergone aprearranged form of convolutive mixing can be used to produce H(z), andthen H(z), the paraconjugate of H(z), can be compared with theprearranged form to see how faithfully the former reconstructs thelatter. Moreover, it can be shown that under certain circumstances H(z)can be used to determine further information about original signals,such as arrival direction and frequency response.

Referring now to FIG. 3, a flow diagram is shown illustrating aniterative process of the invention to produce a strongly decorrelatingpolynomial paraunitary matrix H(z) as the product of a series ofelementary paraunitary matrices. At the beginning of each successiveiteration of the process there is a respective current value h(z) of theevolving paraunitary matrix, i.e. the product of all elementaryparaunitary matrices calculated up to the present time. There is also aset of current signals derived using the most recently calculatedelementary paraunitary matrix incorporated in h(z): b(z) represents thetransformation which converts the original input mixed signals to obtainthe current signals: it is improved repeatedly until applying a furtherelementary paraunitary matrix results in no significant improvement inDCF, at which point h(z) has become H(z).

The process starts at 60, where the initial current signals are theoriginal input signals, and the initial current value of the paraunitarymatrix h(z) is the identity matrix I. At 61 the current signals are usedto produce a current elementary paraunitary matrix; at 62 this matrix isapplied to the current signals to produce new signals. At 63 the newsignals are tested to see whether their strong decorrelation propertieshave improved: if such an improvement has been obtained, at 64 thecurrent elementary paraunitary matrix is used to pre-multiply h(z) toprovide the latest current value of the evolving paraunitary matrix andthe new signals become the latest current signals. The latest currentvalue of h(z) and the latest current signals are then fed back to stage61 for the procedure of stages 61, 62 and 63 to be repeated.

If at 63 the new signals show no improvement over the current signals,the process is considered to have converged to a solution and itterminates at 65: the current value of h(z) is output as H(z) and thecurrent signals are output as the required strongly decorrelatedsignals. Thus the use of elementary paraunitary matrices divides theproblem of strongly decorrelating signals and obtaining H(z) into anumber of sub-problems.

Methods for finding each elementary paraunitary matrix at 61 and fordeciding at 63 whether or not the new signals are improvements over thecurrent signals will now be described. The second relies upon a measureof strong decorrelation of signals being used. An accurate andcomprehensive decorrelation cost function (DCF) is impossible tocalculate from sampled signals, but it is possible to calculate variousdifferent measures of strong decorrelation. Some examples of gooddecorrelation cost functions will now be given. They are based upon thecorrelation terms r(x₁, x₂) of a set of signals defined by:r(x ₁ ,x ₂)=E[x ₁ x ₂]  (18)where E is the expectation operator and x₁ and x₂ are signals in theset: they may be different signals, or the same signals, and they mayhave delays applied to them.

Decorrelation cost functions, denoted by N₁ ^(L), are defined as the sumsquared of all possible correlation terms in a set of signals. In thepresent example correlations are derived between each signal in the setand a relatively delayed replica of itself, i.e. there are nocorrelations between different signals: x₁ and x₂ in Equation (18) aretherefore the same signal with a relative delay therebetween, andcorrelations are obtained for each of a set of delays up to a maximumdelay of L. The maximum delay L can be regarded as an indication of thelocality of the cost function. If L=0 the cost function is referred toas a pointwise cost function, if L is small, i.e. L<5, the cost functionis referred to as a local cost function, otherwise it is a wide costfunction.

Changing the L value changes the properties of N₁ ^(L), but for allvalues of L this is a sensible decorrelation cost function. Inparticular, the two cases of L=0 and L=max (i.e. the largest value thatcan be measured using the sampled signals) have been successfully used.The L=0 DCF leads to an algorithm which imposes spectral majorisation,while the L=max DCF leads to an algorithm that only imposes wide sensedecorrelation, and not spectral majorisation. The process of theinvention may use any convenient DCF, provided it is an accurate DCF andis used consistently throughout: specific problems may suggest specificDCFs.

At 63 current signals and new signals are compared by calculating theirDCFs: the signals having the higher DCF are accepted as being the betterof the two.

To obtain the latest elementary paraunitary matrix at 61, it isnecessary to derive two parameters which characterise such a matrix, adelay parameter d and a rotation parameter θ. Thus identification of theelementary paraunitary matrix reduces to finding (d,θ). The parameter dis a vector parameter, and θ can be a vector parameter if there are morethan two signals. However in this example only the case of two signalsis being considered, thus θ is not a vector.

To find the delay parameter d, a set D of possible values for it isselected. A suitable set for the two-signal case is shown below inEquation (22), i.e. D may be represented as: $\begin{matrix}{D = \left( {\begin{pmatrix}0 \\A\end{pmatrix},\begin{pmatrix}0 \\{A - 1}\end{pmatrix},{\ldots\quad\begin{pmatrix}0 \\1\end{pmatrix}},\begin{pmatrix}0 \\0\end{pmatrix},\begin{pmatrix}1 \\0\end{pmatrix},{\ldots\quad\begin{pmatrix}{A - 1} \\0\end{pmatrix}},\begin{pmatrix}A \\0\end{pmatrix}} \right)} & (19)\end{matrix}$where dots indicate terms not shown, each inner pair of parenthesesrepresents a delay vector δ_(i) having two elements, 0 in an upper orlower position in parentheses indicates an undelayed upper or lowersignal channel respectively, and a non-zero value in such a positionindicates a delayed channel located likewise: each non-zero valuerepresents delaying the associated channel in each case relative to theother channel by an integer value not greater than A. Delays are inunits of a time interval τ between successive signal samples, e.g. adelay (A−1) implies (A−1)τ. There are also other possibilities for D,but before major departures from this scheme are implemented they shouldbe assessed for validity. The number of elements or delay values in D isS, i.e. |D|=S.

Referring now to FIG. 4, there is shown a flow diagram indicating howelementary paraunitary matrices can be found. For all delay vectors δ₁to δ_(S) in D, respective elementary paraunitary matrices are derived,and then the DCFs of the signals these matrices create are calculated.The elementary paraunitary matrix corresponding to the best DCF is thenselected for use. In the drawing, in the two-channel case a pair ofcurrent signals are input at 70. At 71 the signals are replicated Stimes, i.e. as many times as there are delay values in D. In the twochannel case each replica is a pair of current signals: each replicapair is input to a respective signal channel 72 _(i) of which there areS, i.e. i=1 to S. Three channels are shown, 72 ₁, 72₂ and 72 _(S),others being implied by chain lines such as 73. At 74 ₁ to 74 _(S), eachof the channels 72 ₁ to 72 _(S) applies a respective delay vector δ₁, .. . . δ_(S) from D to its replica, i.e. it delays one of the pair ofcurrent signals in its replica relative to the other: the ith channel 72_(i) applies delay vector δ_(i) to the ith replica, and i=1 to S. Ineach of the channels 72 ₁ to 72 _(S) the relative delay produces a pairof signals that are different from those input at 70, although they areno more strongly decorrelated than before.

To make these signals as decorrelated as possible in a pointwise senseusing only a rotation can be achieved by using an instantaneousalgorithm. A review of many methods of imposing instantaneousdecorrelation, both with and without using energy preserving techniques,appears in “Matrix Computations”, John Hopkins University Press 1989 byG. H. Golub and C. F. Van Loan. This covers a wide range of techniquesdeveloped over the last 150 years, and discusses numerical issues causedby the introduction of computers into the field.

Suitable rotations may be calculated by one of several differenttechniques. For many choices of DCF it may be possible to calculate therotation that explicitly maximises the DFC of the signals afer theapplication of the rotation. It may also be possible to try to find thisrotation by a numerical approximation method. A third possibility is touse a pointwise decorrelation algorithm to find the rotation thatimposes pointwise decorrelation; this should lead to an improvement inDCF to at least a predominant part of the maximum improvement attainablein this regard. These three rotation finding possibilities may coincide,as they do for the N₁ ^(L) measure with L=0.

Each channel 72 _(i) now has parameters associated with it consisting ofa delay vector δ_(i) and a rotation angle θ_(i), just as if anelementary paraunitary matrix implementing those parameters had beenapplied to the signals input at 70: moreover, when applied to inputsignals this matrix provides a rotation which produces output signalsthat are as decorrelated as they can be for the relevant value of delayelement δ_(i). The first two elements 74 _(i) and 75 _(i) of eachchannel 72 _(i) therefore in combination simulate an elementaryparaunitary matrix with a single set of delays and a respectiveparameter set (d,θ)=(δ_(i),θ_(i)).

The next stage is to determine which of the channels 72 ₁ to 72 _(S)simulates an elementary paraunitary matrix producing the highestimprovement in strong decorrelation. This is carried out at 76 ₁ to 76_(S), where DCFs of signals output from stages 75 ₁ to 75 _(S)respectively are calculated. These DCFs are compared with one another at77, and whichever channel produces the greatest DCF provides the highestimprovement in strong decorrelation and simulates the best elementaryparaunitary matrix: this matrix is selected as the output of the processindicated in FIG. 4 for use at 61 in FIG. 3.

On making certain assumptions regarding the input signals it is possibleto avoid most of the calculations in this stage. It may be possible topredict the DCFs of signal outputs from stages 75 ₁ to 75 _(S). If thisis possible then it is not necessary to carry out the calculations forany of the channels except the one that will lead to the greatestincrease in DCF. In general if this approach is valid, it is worth doingas the computational savings will be considerable. An example of thistechnique is now given. If the cost function is N₁ ^(L) with L=0 thenthe maximum possible gain in DCF achievable for a given delay is simplythe correlation between the two signals after that delay has beenapplied. By using the circulant assumption, this correlation can becalculated for all possible delays by using Fourier transforms. This isa prior art technique for calculating the cross correlation between twosignals, and is computationally inexpensive compared with calculatingthe cross correlations for each delay individually. Thus the delay whichleads to the maximum cross correlation between the two signals can beselected, without having to apply delays and calculate rotations foreach of the channels.

It is not in fact essential to choose whichever channel produces thegreatest DCF: another channel giving a value relatively close to thisvalue would be acceptable. Also a channel might be selected to give amore convenient value for one of the parameters, e.g. a shorter delay.If however a channel is chosen giving a value insufficiently close tothe greatest DCF, then the process of signal separation may fail becauseit might be impossible to rectify matters in subsequent iterations. So achannel should be selected which gives a DCF or improvement indecorrelation to at least a predominant part of the maximum in thisregard obtainable over the ranges of delay and rotation parametersemployed.

The foregoing embodiment of the invention related to imposing wide sensedecorrelation in a two-channel case, i.e. two signal sources, tworeceivers and two paths from each source to receivers. If there are morethan two channels further processing is necessary. However the procedureof the invention remains fundamentally the same: i.e. find a paraunitarymatrix to impose strong decorrelation on signals in the channels, thematrix being a product of successive elementary paraunitary matriceseach of which maximises a measure of strong decorrelation of outputsresulting from its application to input signals.

There are several ways of extending the invention to more than twochannels. One is a simple, direct extension of the previous embodiment.Others use signal pairs, attempt to impose pairwise independence, andsweep through all signal pairs in some order.

FIG. 5 shows the decomposition of a paraunitary matrix with degree N forthe multi-channel case (i.e. three or more channels) as given by theVaidyanathan prior art previously mentioned. It is very similar to thetwo-channel case of FIG. 1. The difference is that Givens rotations,which are essentially two-channel operators, have been replaced withmore general rotations operating on all channels at once. This can beexpressed mathematically as:—H _(N)(z)=R _(N) . . . Λ(z)R ₁Λ(z)R ₀  (20)where as before Λ(z) is a matrix representing successive unit delays inthe lowermost channel, and R_(i) (i=0 to N) represents the ith rotation.In the n channel case each rotation has n(n−1)/2 parameters, whichreduces to one parameter when n=2 for Givens rotations used in thetwo-channel case.

The parameterisation of Equation (20) involves the same problems as inthe two-channel case described with reference to FIG. 1. Here again theprocedure used to deal with the problems is the decomposition of aparaunitary matrix into elementary paraunitary matrices: the elementaryparaunitary matrices are found by considering the decomposition andconstraining all but one rotation. One rotation is unconstrained, butall others are constrained to either leave inputs unchanged or toreorder channels so a different channel receives a subsequent delay.Thus all rotations apart from one are constrained to be permutations. Inthe present example of an elementary paraunitary matrix with only oneset of delays, the unconstrained rotation is the last rotation, and thematrix can thus be expressed as:—V _(d,θ)(z)=R _(θ)Λ(z) . . . P ₁Λ(z)P ₀ =R _(θ) D _(d)(z)  (21)

Here P_(i) is the ith permutation matrix and R_(θ) is the rotationmatrix parameterised by a vector of n(n−1)/2 rotation angles θ. Thecombination of all permutation matrices and delay matrices is denoted byD_(d), which is referred to as the delay matrix. Its index, d, is alength n vector whose n elements are delays applied in respectivechannels.

Equation (21) shows the extended form of an elementary paraunitarymatrix. First a delay matrix is applied, which consists of differingdelays being applied in all the channels. Secondly a rotation isapplied. As before this is a good choice for an elementary buildingblock, as it contains only one operation that mixes the signals.However, this form is parameterisable, the parameters being the vectorsd and θ. This means that the decomposition of a paraunitary matrix intoelementary paraunitary matrices can be written as:—H(z)=V _(d) _(N) _(,θ) _(N) (z) . . . V _(d) ₁ _(,θ) ₁ (z)V _(d) ₀ _(,θ)₀ (z)  (22)

The methodology for finding this expansion is largely the same as in thetwo-channel case. It only requires modification in selecting D, i.e. thedelay set or set of possible delays for forming the first of the seriesof elementary paraunitary matrices. In the two-channel case the elementsof D were treated as delays applied to only one channel. In the moregeneral case delays are allowed in more than one channel. As in thetwo-channel case there is a countable infinity of possibilities, so asensible finite set has to be selected to form D.

The number of elements in D could be fixed by setting an upper limitdenoted by l to the number of delays in it. Thus D can consist of allpossible ways of allocating up to l delays between n channels. Applyingthe same delay to all channels at the same time is equivalent toapplying no delay at all, so all elements in D should have at least onechannel that is not delayed. The number of elements in D can grow veryfast as n and l increase. As before this scheme for D is merelyselection of a sensible option. Other possibilities exist but usuallyneed more justification to ensure validity. It is possible to justifysome choices (such as allowing delays to be applied to only one channel)by saying they allow a higher l for a fixed size of D. Thus they allowthe algorithm to consider longer delays, at the cost of not allowing itto implement more complicated combinations of smaller delays.

However it is defined, once D is fixed the process of the invention canproceed to impose wide sense decorrelation in the multi-channel case ina similar way to the two-channel case of the previous embodiment. Thereare two differences: firstly, the two-element delay vectors δ_(i) inFIG. 4 are replaced by delay vectors (elements of the delay set D)having n elements, where n is the number of channels. This allows thealgorithm to apply an element of D to the signals. Secondly, therotations calculated by known pointwise decorrelation algorithms atstages 75 ₁ to 75 _(S) are now n channel rotations. These can beobtained by successively applying n(n−1)/2 two channel rotations, and socan be parameterised by n(n−1)/2 angle parameters: Prior art pointwisedecorrelation algorithms previously mentioned provide current signalsand a rotation matrix which generates them, so the process will not bedescribed further.

This approach to the multi-channel case may give difficulty over thesize of the delay set D: because each element in D is a vector, for eachelement a whole n by n pointwise decorrelation algorithm has to becarried out. As the delay set upper limit l grows this becomes thedominant time constraint and slows processing down. So-called “sweeping”algorithms were developed in the prior art to attempt to mitigatesimilar problems in the instantaneous case.

The principle of using energy preserving transforms to zero thecorrelation between a set of signals has a long history. The firstreference is “Uber ein Leichtes Verfahren Die in der Theorie derSacularstroungen Vorkommendern Gleichungen Numerisch Aufzulosen”Crelle's J. 30 p 51-94, 1846 by C. G. J. Jacobi. He introduced the ideaof tackling the multichannel case by sweeping through pairs of signalsin such a way as to strictly decrease the total measure of correlationbetween different signals in the problem at each stage. This techniqueis known as using Jacobi sweeps

It has been shown that cycling through different signal pairings andimposing pairwise decorrelation leads to good decorrelation of severalsignals, provided that there are repeated sweeps through all availablesignal pairs. The technique for decorrelating several signals byrepeated sweeps through all available signal pairs imposing pairwisedecorrelation on each pair is known as the Jacobi algorithm. At any onestage it only works on making two signals pairwise decorrelated.

This example of the invention provides a procedure called the sweepingprocess for finding an elementary paraunitary matrix based algorithmthat follows the same approach used in the Jacobi algorithm. Theprocedure has a prearranged ordering of all available signal pairingsthrough which it works in a sweep. Each signal pair is processed usingthe two-channel process of the first embodiment of the invention. Theremaining signals are aligned in time (delayed) so they have at leastapproximately the same delay applied to them as the signal pair beingprocessed. At the end of a sweep the same termination condition from theprevious embodiment of the invention can be applied to decide if anothersweep needs to be carried out or if no more significant improvement isobtainable.

FIG. 6 illustrates this procedure for imposing wide sense decorrelationon three signals A, B and C. As mentioned earlier the process of theinvention implements the Sequential Best Rotation or SBR algorithm. Inthe drawing, “SBR2” indicates application of the SBR algorithm for atwo-channel case as in the previous embodiment.

In an initial phase, signals A and B are processed with SBR2 at 90, andat 91 signal C is aligned in time for substantial synchronisation withthe output from 90. In this example the elementary paraunitary matricesimplement only one set of delays, and time alignment at 91 is carriedout by delaying C by half the order of the polynomial matrix determinedat 90. In other words C is delayed in time by half the sum total of thedelays applied in the processing of A and B irrespective of whichchannel the or (as the case may be) each delay lies in. Alignment of Cin this way is not essential but has been found to give acceptableresults. It is not necessary to allow for processing delay to implementSBR2 at 90, because the signals A, B and C are sequences of digital datasamples stored when not in use: “delaying” C relative to A or B merelymeans changing the sample number in a sequence input to a subsequentprocessing stage.

At stages 90 and 91 in combination, at least one elementary paraunitarymatrix is applied to all three signals: such a matrix can be generatedby extending the or as the case may be each two channel elementaryparaunitary matrix produced at 90 to cover three channels. The or eachmatrix applies rotation and delay to signals A and B but it applies adelay to signal C without rotating it. If more than one elementaryparaunitary matrix is determined at 90, all delays applied to signal Cmust add up to the total delay applied to C at 91.

In a second phase, signal B (after processing at 90) and signal C (afterdelay at 91) are processed with SBR2 at 92, and signal A (afterprocessing at 90) is aligned in time at 93 with the output from 92. In athird phase, signal A (after alignment at 93) and signal C (afterprocessing at 92) are processed with SBR2 at 94, and signal B (afterprocessing at 92) is aligned in time at 95 with the output from 94.

The next stage at 96 is to determine whether or not the signals havebecome sufficiently strongly decorrelated to warrant terminating theprocess, or whether it is necessary to continue: as before thistermination condition is satisfied if DCFs of the signals have notimproved significantly: if so, processing terminates; if not arespective elementary paraunitary matrix determined at each of 90, 92and 94 is used to pre-multiply a stored product of those preceding itand the latest resulting product replaces that previously stored; thelatest current signals from SBR2 94 and alignment 95 are then fed backto 90 and 91 via a feedback loop 97 for the procedure to be repeated.

In this embodiment, one or more two-channel elementary paraunitarymatrices may be applied in the SBR2 stages at 90, 92 and 94 for eachcirculation round the loop of items 90 to 97. It is possible to apply asmany elementary paraunitary matrices at these stages as one wishes, upto the number necessary to meet the termination condition for the SBR2algorithm. It is preferred to apply only one elementary paraunitarymatrix at the stages 90, 92 and 94 because it reduces the need forsignal alignment at 91, 93 and 95: it uses the philosophy of applyingthe more important rotations first, because earlier rotations in theSBR2 process tend to give more performance gain than later equivalents.

If, rather than operating on signal pairs in a fixed order, the pairthat leads to the greatest DCF increase is operated on first, anotherapproach for imposing wide sense decorrelation can be produced. Thisapproach also operates on pairs, but pair ordering is not fixed. Insteadall possible pairs are considered and the one that leads to the greatestimprovement in DCF is accepted. This is known as a searching process, asit searches through possible signal pairs. FIG. 7 illustrates theprocedure for three signals A, B and C. The signals are replicatedn(n−1)/2 times where n is the number of signals: n=3 in the presentexample, which coincidentally means there are three replicas of thesignals. In an initial phase, using a first replica of the three signalsA, B and C, at 100 signals A and B are processed with SBR2, and, in thisexample using elementary paraunitary matrices with one set of delays,signal C is aligned in time with the output so produced. Here again a(or as the case may be each) two channel elementary paraunitary matrixproduced by the application of SBR2 to two signal channels is extendedto an n channel elementary paraunitary matrix that delays (but does notrotate) the or each remaining signal C to preserve time alignment withSBR2-processed signals A and B.

For a second replica, at 101 signals A and C are processed with SBR2,and signal B is aligned with the output. For a third replica, at 102signals B and C are processed with SBR2 at 94, and signal A is alignedwith the output. At 103, the DCFs of the three output signal groups aredetermined: the group with the greatest DCF is selected, and if the DCFshows an improvement over the input signals to this stage, then theoperation is accepted.

As before, if the DCF of all the signals has been improved, processingcontinues with new current signals and new current matrix being fed backvia a loop 104 to 100, 101 and 102. If there has been no significantimprovement, processing terminates and at 105 the current signals andparaunitary matrix are output.

Again the SBR2 steps 100, 101 and 102 may consist of applying a fullSBR2 process of calculating a series of elementary paraunitary matricesuntil no improvement is obtained, or just finding one elementaryparaunitary matrix. The latter case is equivalent to a restrictedversion of a full SBR: i.e. ignoring time alignment, the delays arerestricted to only delaying one channel and the rotations to onlyapplying a rotation to the delayed channel and one other. This reducesthe size of the delay set and hence increases the speed of computationof elementary paraunitary matrices.

Although the elementary paraunitary matrices described with reference toFIG. 2 are sufficient for implementing the invention, practicalapplications of the invention may work better if an extended definitionof elementary paraunitary matrices is used. A more general definition ofan elementary paraunitary matrix is a polynomial matrix which representsapplying an initial set of delays to input signals, then applying asingle unitary transformation (rotation) to the delayed input signals,followed by applying a further set of delays to the delayed and rotatedinput signals. The initial set of delays is defined as the leadingdelays and the further set of delays as the terminal delays.

FIG. 8 illustrates a single elementary paraunitary matrix implementation120 of this new form: as in FIGS. 2 and 3 this implementation is shownfor a two channel case for the sake of simplicity. In the implementation120, signals in upper and lower channels 121 and 122 pass into a matrixof four ganged single pole double throw switches 123 controlling routingor otherwise via a delay cell z^(di). The switches 123 are all either inthe UP position or in the DOWN position: when DOWN, upper channelsignals are delayed at Z^(di) but lower channel signals are not, andwhen UP the reverse is the case. This enables either channel 121 or 122to be delayed relative to the other.

On leaving the switches 123, an upper channel signal is multiplied byGivens rotation parameters −s_(i) and c_(i) at amplifiers 124 and 125respectively. Similarly, on leaving the switches 123, a lower channelsignal is multiplied by Givens rotation parameters s_(i) and c_(i) atamplifiers 126 and 127 respectively. Here s_(i) and c_(i) arerespectively the sine and cosine of a rotation angle θ_(i) implementedat 120. Each c_(i) product is summed with the s_(i) or −s_(i) productinvolving the signal from the other channel: this provides sums whichpass along respective channels 128 and 129 to another matrix of fourganged single pole double throw switches 130 controlling routing orotherwise via a delay cell Z^(Δi). These operate in a similar way to theswitches 123.

A paraunitary matrix can be constructed by replacing each of the dottedline blocks V₀(z) to V_(N)(z) in FIG. 2 by a respective implementation120 of the new elementary paraunitary matrix: each of these blocks 120or V_(i)(z) (i=0, 1, to N) represents a respective single elementaryparaunitary matrix that implements a respective leading delay d_(i),Givens rotation through an angle θ_(i) and terminal delay Δ_(i).

A full paraunitary matrix parameterisation in accordance with theinvention in terms of the new elementary paraunitary matrices can beexpressed as follows:—H _(N)(z)=V _(d) _(N) _(,θ) _(N) _(,Δ) _(N) (z) . . . V _(d) ₁ _(,θ) ₁_(,Δ) ₀ (z)V _(d) ₀ _(,θ) ₀ _(,Δ) ₀ (z)  (23)

Although the full specification of an elementary paraunitary matrixincludes an arbitrary terminal set of delays parameterised by Δ_(i), theinvention works best if one of two specific types of elementaryparaunitary matrix are used. Type 1 elementary paraunitary matrices,which were described with reference to FIG. 2, have no (or zero)terminal delays, i.e. Δ_(i)=0. Type 2 elementary paraunitary matriceshave terminal delays which apply an inverse of respective leadingdelays, i.e. Δ_(i)=−d_(i). It will be recalled that a negative delay inone channel is implemented by a positive delay in the or each (as thecase may be) other channel.

The process of the invention using the new or type 2 elementaryparaunitary matrices, to find an appropriate paraunitary matrix thatimposes wide sense decorrelation on the signals is very similar to thatpreviously described. The appropriate paraunitary matrix is notparameterised as disclosed by Vaidyanathan, but it is insteadparameterised in terms of elementary paraunitary matrices (see Equation(23)) calculated sequentially. Once derived, each elementary paraunitarymatrix is applied to the signals used to calculate it: this transformsthe signals into a form suitable for calculating the next such matrix.

A delay pair is defined as a pair of delays, one the leading delay andthe other the terminal delay, used in an elementary paraunitary matrix.To calculate an elementary paraunitary matrix a set of possible delaypairs is considered. A latest elementary paraunitary matrix is chosenwhich maximises the DCF calculated for signals transformed byapplication of all such matrices calculated up and including the latest.To identify the latest elementary paraunitary matrix its delay androtation parameters are required. There is a countable infinity ofdifferent possible delays, but an uncountable number of differentrotation parameters. Realistically the countable infinity of possibledelays will be reduced to a finite number of different delays to beconsidered, either by constraining them or by making the assumption ofcirculant statistics and hence constraining them to be less than thedata length (number of samples of each signal used in decorrelation). Iftype 1 or type 2 elementary paraunitary matrices (as defined above) areto be used the terminal delays are fixed by the value of the leadingdelays.

Once a trial value of leading delay has been applied to the signals,rotation parameters s_(i) and c_(i) are chosen. These parameters couldbe chosen by using a technique which explicitly maximises the DCF ofoutput signals.

It is therefore possible to look at a set of trial delay pairs and foreach delay pair to attach the rotation that maximised the DCF after theapplication of the leading delay, the rotation and the terminal delay.This produces a set of possible elementary paraunitary matrices, eachassociated with a respective DCF found by applying it to signals andfinding the DCF of the resulting output. Under certain assumptions itmay be possible to avoid applying the elementary paraunitary matrices tothe signals, and instead to calculate the DCF of the resulting outputdirectly. This has computational advantages, but is only possible insome cases.

The elementary paraunitary matrix that is selected is that having thebest DCF output. This stage is repeated until there is no gain in DCF tobe found, or the gain in DCF is regarded as being too small to warrantthe extra complexity in the strongly decorrelating matrix.

The iterative process for calculating a strongly decorrelatingparaunitary matrix H(z) as shown in FIG. 3 remains nearly the same forthe new or type 2 elementary paraunitary matrices. The only differenceis that to obtain the latest elementary paraunitary matrix at 61, it isnecessary to derive parameters which characterise such a matrix, delayparameters d and Δ and a rotation parameter θ. Thus identification ofthe elementary paraunitary matrix reduces to finding (d,θ,Δ). Theparameters d and Δ are vector parameters, and θ can be a vectorparameter if there are more than two signals. In the case of type 2elementary paraunitary matrices, Δis entirely dependent upon d, and soonce d has been found (as described earlier) Δ has also been found.

The additional parameter Δ results in a small change as compared to theprocess described with reference to FIG. 4. This change is illustratedin FIG. 9, which shows a processing channel 140 _(i) (i=1 to S) toreplace channel 72 _(i) in FIG. 4. Ignoring the fact the delay at 74_(i) has been called the leading delay at 142 _(i), the only differencebetween the channels 72 _(i) and 142 _(i) is that the latter includesprovision for application of a terminal delay Δ_(i) at 144 _(i) afterrotation at 146 _(i) and before DCF calculation at 148 _(i).

The equations and procedures given in the foregoing description canclearly be implemented by an appropriate computer program comprisingprogram instructions on an appropriate carrier medium and running on aconventional computer system. The carrier medium may be a memory, afloppy or compact or optical disc or other hardware recordal medium, oran electrical signal. Such a program is straightforward for a skilledprogrammer to implement from the foregoing description without requiringinvention, because it involves well known computational procedures

1. The present invention provides a method of strong decorrelation ofsignals which are input, the method including the steps of:— a)processing the input signals to determine delay and rotation parameterswhich implement at least one elementary paraunitary matrix and transformthe input signals into output signals to obtain improvement in a measureof strong decorrelation; b) assessing the improvement in the measure ofstrong decorrelation, and c) if the improvement is significantdesignating the output signals as input signals and iterating steps a)and b), and if the improvement is not significant designating the outputsignals as signals decorrelated in a wide sense.
 2. A method accordingto claim 1 wherein the delay and rotation parameters which transform theinput signals characterise a single elementary paraunitary matrix.
 3. Amethod according to claim 2 including producing a paraunitary matrix bycumulatively multiplying successive elementary paraunitary matricesproduced by iterating step a).
 4. A method according to claim 2 whereinthe range of signal delay parameters is a set of discrete delay vectors,and the delay and rotation parameters are determined by generating arespective version of the input signals delayed by each delay vector inthe set, and for each version finding rotation parameters which at leastapproach producing maximisation of output signals' strong decorrelation.5. A method according to claim 4 wherein the rotation parameters whichat least approach producing maximisation of output signals' strongdecorrelation are determined using an algorithm for pointwisedecorrelation of the kind used in instantaneous decorrelation.
 6. Amethod according to claim 1 involving n input signals where n is aninteger greater than 2, wherein the range of signal delay parameters isa set of n-element delay vectors and the range of signal rotationparameters is a set of n(n−1)/2 angle parameters.
 7. A method accordingto claim 1 involving n input signals where n is an integer greater than2, wherein step a) comprises determining delay and rotation parameterswhich implement at least one elementary paraunitary matrix providing forrotation of a pair of input signals and relative delay of the or as thecase may be each other input signal.
 8. A method according to claim 7wherein the n input signals are associated with respective channelswherein step a) has n(n−1)/2 successive stages each associated with atleast one respective elementary paraunitary matrix and each providingfor rotation of signals associated with a respective pair of channelsand provision of relative delay associated with the or as the case maybe each other channel, the first stage is arranged to process the inputsignals and the or as the case may be each subsequent stage is arrangedto receive signals processed in the respective preceding stage.
 9. Amethod according to claim 1 involving a set of n input signals where nis an integer greater than 2, the method having the steps of: a)producing n(n−1)/2 replicas of the set of input signals, b) in eachreplica selecting a respective signal pair differing to that selected inother replicas, and c) the step of processing the input signals todetermine delay and rotation parameters being carried out for eachreplica and comprising: i) determining delay and rotation parameterswhich implement at least one elementary paraunitary matrix providing forrotation of the respective selected signal pair only, and ii)determining which replica when transformed by the associated at leastone elementary paraunitary matrix gives rise to transformed signalscorresponding to improvement in a measure of decorrelation by at least amajor part of a maximum extent obtainable over the replicas anddesignating these transformed signals as output signals.
 10. A methodaccording to claim 1 wherein the at least one elementary paraunitarymatrix implements at least one leading delay, rotation and terminaldelay
 11. Computer apparatus for strong decorrelation of signals, theapparatus being programmed for reception of input signals and also beingprogrammed:— a) to process the input signals to determine delay androtation parameters which implement at least one elementary paraunitarymatrix and transform the input signals into output signals to obtainimprovement in a measure of strong decorrelation; b) to assess theimprovement in the measure of strong decorrelation, and if theimprovement is significant to designate the output signals as inputsignals and iterate a) and b), and if the improvement is not significantto designate the output signals as signals decorrelated in a wide sense.12. Apparatus according to claim 11 wherein the delay and rotationparameters which transform the input signals characterise a singleelementary paraunitary matrix.
 13. Apparatus according to claim 12characterised in that the computer equipment is programmed to produce aparaunitary matrix by cumulatively multiplying successive elementaryparaunitary matrices produced in iterative processing.
 14. Apparatusaccording to claim 12 wherein the range of signal delay parameters is aset of discrete delay vectors, and the computer apparatus is programmedto determine the delay and rotation parameters by generating arespective version of the input signals delayed by each delay vector inthe set, and to find for each version rotation parameters which at leastapproach producing maximisation of output signals' strong decorrelation.15. Apparatus according to claim 14 programmed to determine the rotationparameters which at least approach producing maximisation of outputsignals' strong decorrelation using a pointwise decorrelation algorithm.16. Apparatus according to claim 11 programmed to receive n inputsignals where n is an integer greater than 2, also programmed todetermine delay and rotation parameters which implement at least oneelementary paraunitary matrix providing for rotation of a pair of inputsignals and relative delay of the or as the case may be each other inputsignal.
 17. Apparatus according to claim 16 programmed to definerespective channels for the n input signals to process the input signalsin n(n−1)/2 successive stages each associated with at least onerespective elementary paraunitary matrix and each providing for rotationof signals associated with a respective pair of channels and provisionof relative delay associated with the or as the case may be each otherchannel, the first such stage involving processing the input signals andthe or as the case may be each subsequent stage involving processingsignals resulting from processing in the respective preceding stage. 18.Apparatus according to claim 11 programmed to receive a set of n inputsignals where n is an integer greater than 2, and also programmed to: a)produce n(n−1)/2 replicas of the set of input signals, b) in eachreplica select a respective signal pair differing to that selected inother replicas, and c) implement processing of the input signals todetermine delay and rotation parameters for each replica as inputsignals and determine: i) delay and rotation parameters which implementat least one elementary paraunitary matrix providing for rotation of therespective selected signal pair only, and ii) which replica whentransformed by the associated at least one elementary paraunitary matrixgives rise to transformed signals corresponding to improvement in ameasure of strong decorrelation by at least a major part of a maximumextent obtainable over the replicas and designate these transformedsignals as output signals.
 19. Apparatus according to claim 11 whereinthe at least one elementary paraunitary matrix implements at least oneleading delay, rotation and terminal delay.
 20. A computer programmehaving instructions for implementing strong decorrelation of signalsinput to computer apparatus, the computer programme having instructionsfor controlling the computer apparatus:— a) to process the input signalsto determine delay and rotation parameters which implement at least oneelementary paraunitary matrix and transform the input signals intooutput signals to obtain improvement in a measure of strongdecorrelation; b) to assess the improvement in the measure of strongdecorrelation, and if it is significant to designate the output signalsas input signals and iterate a) and b), and if the improvement is notsignificant to designate the output signals as signals decorrelated in awide sense.
 21. A computer programme according to claim 20 wherein thedelay and rotation parameters which transform the input signalscharacterise a single elementary paraunitary matrix.
 22. A computerprogramme according to claim 21 having instructions for controllingcomputer apparatus to implement the step of producing a paraunitarymatrix by cumulatively multiplying successive elementary paraunitarymatrices produced by iterating processing of the input signals todetermine delay and rotation parameters.
 23. A computer programmeaccording to claim 21 wherein the range of signal delay parameters is aset of discrete delay vectors, and the computer programme havinginstructions for controlling computer apparatus to provide for the delayand rotation parameters to be determined by generating a respectiveversion of the input signals delayed by each delay vector in the set,and for each version finding rotation parameters which at least approachproducing maximisation of output signals' strong decorrelation.
 24. Acomputer programme according to claim 23 having instructions forcontrolling computer apparatus to provide for the rotation parameterswhich at least approach producing maximisation of output signals' strongdecorrelation to be determined using a pointwise decorrelationalgorithm.
 25. A computer programme according to claim 20 havinginstructions for controlling computer apparatus to receive n inputsignals where n is an integer greater than 2, and to provide forprocessing the input signals to determine delay and rotation parametersto comprise determining such parameters which implement at least oneelementary paraunitary matrix providing for rotation of a pair of inputsignals and relative delay of the or as the case may be each other inputsignal.
 26. A computer programme according to claim 25 havinginstructions for controlling computer apparatus to define respectivechannels for the n input signals and also having instructions forcontrolling computer apparatus to process the input signals to determinedelay and rotation parameters to have n(n−1)/2 successive stages eachassociated with at least one respective elementary paraunitary matrixand each providing for rotation of signals associated with a respectivepair of channels and provision of relative delay associated with the oras the case may be each other channel, the first stage being arranged toprocess the input signals and the or as the case may be each subsequentstage being arranged to receive signals processed in the respectivepreceding stage.
 27. A computer programme according to claim 22 havinginstructions for controlling computer apparatus to receive a set of ninput signals where n is an integer greater than 2, and furtherinstructions for controlling the computer apparatus to: a) producen(n−1)/2 replicas of the set of input signals, b) in each replica selecta respective signal pair differing to that selected in other replicas,and c) carry out processing of determine delay and rotation parametersfor each replica as input signals by: i) determining delay and rotationparameters which implement at least one elementary paraunitary matrixproviding for rotation of the respective selected signal pair only, andii) determining which replica when transformed by the associated atleast one elementary paraunitary matrix gives rise to transformedsignals corresponding to improvement in a measure of strongdecorrelation by at least a major part of a maximum extent obtainableover the replicas and designating these transformed signals as outputsignals.
 28. A computer programme according to claim 22 wherein the atleast one elementary paraunitary matrix implements at least one leadingdelay, rotation and terminal delay.